﻿using System;

namespace LxCFD
{
    public class BiCGStabInZone : BiCGBase
    {
        public double[] r0;
        public double[] r;//s也用r存储
        public double[] p;
        public double[] v;
        public double[] y;
        public double[] z;
        public double[] t;

        public void CreateLocal(IProvideVariableCount variableCount)
        {
            int length = variableCount.GetVariableCount();
            LxTool.CheckCreateArray(ref p, length);
            LxTool.CheckCreateArray(ref v, length);
            LxTool.CheckCreateArray(ref y, length);
            LxTool.CheckCreateArray(ref z, length);
            LxTool.CheckCreateArray(ref t, length);
        }
        public void CreateCommon(double[] rCommon, double[] r0Common)
        {
            this.r = rCommon;
            this.r0 = r0Common;
        }
        public void Initialize0ValueInArray(IndexNumber indexNumber, int taskIndex)
        {
            indexNumber.CellIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                r[i1] = 0;
                r0[i1] = 0;
                p[i1] = 0;
                v[i1] = 0;
                y[i1] = 0;
                z[i1] = 0;
                t[i1] = 0;
            }
        }
        public void Cal_r_Dot_r0(IndexNumber indexNumber, int taskIndex)
        {
            indexNumber.CellIDs.GetSeg(out int i1, out int endp1, taskIndex);
            double temp = 0;
            for (; i1 < endp1; i1++)
            {
                temp += r0[i1] * r[i1];
            }
            ThreadsMath.AddToSum(temp);
        }
        public void Cal_xArray(IndexNumber indexNumber, double[] x, int taskIndex)
        {
            indexNumber.CellIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                x[i1] += alpha * y[i1] + omega * z[i1];
            }
        }
        public void Cal_sArray(IndexNumber indexNumber, int taskIndex)
        {
            Cal_rMinus(indexNumber, BiCGBase.alpha, v, taskIndex);
        }
        public void Cal_rArray(IndexNumber indexNumber, int taskIndex)
        {
            Cal_rMinus(indexNumber, BiCGBase.omega, t, taskIndex);
        }
        protected void Cal_rMinus(IndexNumber indexNumber, double factor, double[] subtrahend, int taskIndex)
        {
            indexNumber.CellIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                r[i1] -= factor * subtrahend[i1];
            }
        }
        public void Cal_pArray(IndexNumber indexNumber, int taskIndex)//2
        {
            indexNumber.CellIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                p[i1] = r[i1] + beta * (p[i1] - omega * v[i1]);
            }
        }
        public void Cal_yArray(IndexNumber indexNumber, double[] ADiag, int taskIndex)
        {
            DInversePrecond(indexNumber, y, ADiag, p, taskIndex);
        }
        public void Cal_zArray(IndexNumber indexNumber, double[] ADiag, int taskIndex)
        {
            DInversePrecond(indexNumber, z, ADiag, r, taskIndex);
        }
        protected void DInversePrecond(IndexNumber indexNumber, double[] left, double[] ADiag, double[] right, int taskIndex)
        {
            indexNumber.CellIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                left[i1] = right[i1] / ADiag[i1];
            }
        }
    }
}

//protected void createMaxLengthArray(ref double[] Array, int length)
//{
//    if (Array == null || Array.Length < length)
//    {
//        Array = new double[length];
//    }
//}
//public void cal_omega(IndexNumber indexNumber, int taskIndex)
//{
//    indexNumber.CellIdx.GetSeg(out int i1, out int endp1, taskIndex);
//    double temp1 = 0, temp2 = 0;
//    for (; i1 < endp1; i1++)
//    {
//        double stemp = r[i1];
//        temp1 += t[i1] * stemp;
//        temp2 += stemp * stemp;
//    }
//    TasksControl.Storage1[taskIndex] += temp1;
//    TasksControl.Storage2[taskIndex] += temp2;
//}
//public static void CalDotProduct(int start, int endp1, double[] array1, double[] array2, int taskIndex)
//{
//    double rhoAtitemp = 0;
//    for (; start < endp1; start++)
//    {
//        rhoAtitemp += array1[start] * array2[start];
//    }
//    TasksControl.Storage1[taskIndex] += rhoAtitemp;
//}
//public static void CalArrayMinus(int start, int endp1, double[] minuend, double factor, double[] subtrahend)
//{
//    for (; start < endp1; start++)
//    {
//        minuend[start] -= factor * subtrahend[start];
//    }
//}
//public static void Calpi(int start, int endp1, double[] p, double[] r, double beta, double omega, double[] v)
//{
//    for (; start < endp1; start++)
//    {
//        p[start] = r[start] + beta * (p[start] - omega * v[start]);
//    }
//}
//public static void Calx(int start, int endp1, double[] xArray, double alpha, double[] p, double omega, double[] s)
//{
//    for (; start < endp1; start++)
//    {
//        xArray[start] += alpha * p[start] + omega * s[start];
//    }
//}
//public void Cal_sArray(IndexNumber indexNumber, int taskIndex)
//{
//    cal_rMinus(indexNumber, alpha, v, taskIndex);
//}
//protected void cal_y(IndexNumber indexNumber, int taskIndex)//3
//{
//    Precond.CalKInverseProduct(indexNumber, y, pAti, taskIndex);
//}
//protected void cal_vAti(IndexNumber indexNumber, int taskIndex)//4
//{
//    AProduct.CalAProduct(indexNumber, vAti, y, taskIndex);
//}
//protected void cal_rAt0ProductvAti(IndexNumber indexNumber, int taskIndex)//5.1
//{
//    indexNumber.CellIdx.GetSeg(out int i1, out int endp1, taskIndex);
//    double temp = 0;
//    for (; i1 < endp1; i1++)
//    {
//        temp += rAt0[i1] * vAti[i1];
//    }
//    TasksControl.Storage1[taskIndex] += temp;
//}
//protected void cal_s(IndexNumber indexNumber, int taskIndex)//6
//{
//    indexNumber.CellIdx.GetSeg(out int i1, out int endp1, taskIndex);
//    for (; i1 < endp1; i1++)
//    {
//        s[i1] = r[i1] - alpha * vAti[i1];
//    }
//}
//protected void cal_z(IndexNumber indexNumber, int taskIndex)//7
//{
//    Precond.CalKInverseProduct(indexNumber, z, s, taskIndex);
//}
//protected void cal_t(IndexNumber indexNumber, int taskIndex)//8
//{
//    AProduct.CalAProduct(indexNumber, t, z, taskIndex);
//}
//protected void cal_KtKsKtKt(IndexNumber indexNumber, int taskIndex)//9.1
//{
//    Precond.CalKInversets(out double fenzi, out double fenmu, indexNumber, t, s, taskIndex);
//    TasksControl.Storage1[taskIndex] += fenzi;
//    TasksControl.Storage2[taskIndex] += fenmu;
//}
//protected void cal_xAti(IndexNumber indexNumber, int taskIndex)//10
//{
//    indexNumber.CellIdx.GetSeg(out int i1, out int endp1, taskIndex);
//    for (; i1 < endp1; i1++)
//    {
//        x[i1] += alpha * y[i1] + omega * z[i1];
//    }
//}
//protected void cal_rAti_1(IndexNumber indexNumber, int taskIndex)//11,但是改为了rAti-1
//{
//    indexNumber.CellIdx.GetSeg(out int i1, out int endp1, taskIndex);
//    for (; i1 < endp1; i1++)
//    {
//        r[i1] = s[i1] - wAti_1 * t[i1];
//    }
//}
//public void Cal_rho(IndexNumber indexNumber, int taskIndex)//1.1
//{
//    cal_DotProduct(indexNumber, r0, r, taskIndex);
//}
//public void Cal_alphaPart1(IndexNumber indexNumber, int taskIndex)
//{
//    cal_DotProduct(indexNumber, r0, v, taskIndex);
//}
//indexNumber.CellIdx.GetSeg(out int i1, out int endp1, taskIndex);
//double rhoAtitemp = 0;
//for (; i1 < endp1; i1++)
//{
//    rhoAtitemp += rAt0[i1] * r[i1];
//}
//TasksControl.Storage1[taskIndex] += rhoAtitemp;